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1.  Introduction 


Early  investigations  on  thermo-mechanical  interaction  resulted  in  the 
classical  studies  of  thermoelasticity.  Two  major  treatises  have  been  published 
which  summarize  the  contributions  in  thermoelasticity  up  to  the  time  of 
early  nineteen  sixties  [1-2] .  Most  commonly,  problems  of  thermal  stress  were 
formulated  such  that  the  temperature  field  does  not  depend  on  the  stress 
field,  whereas  the  stress  field  is  affected  by  the  thermal  expansion  or 
contraction  of  the  material.  However,  the  fully  coupled  system  of  equations 
of  thermoelasticity  containing  coupling  terms  in  both  the  equation  of  motion 
and  the  energy  balance  equation  has  received  relatively  little  effort  until 
more  recently.  This  development  is  largely  due  to  the  fact  thfct  these  equations 
by  and  large  are  not  easily  accessible  by  the  available  analytical  techniques. 
Some  successes  were  achieved  in  isolated  cases  [3-51.  In  these  cases  several 
methods  of  analysis  including  Laplace  transform  or  perturbation  series  have 
been  employed.  For  the  fully  coupled  equations  with  arbitrary  coupling 
coefficient  and  with  general  type  of  boundary  conditions  one  must  resort 
to  numerical  techniques.  Numerical  approach  to  the  solution  of  the  fully 
coupled  thermoelasticity  equation  has  appeared'  in  the  literature  [6] .  More 
recently,  finite  element  method  has  been  applied  in  solving  the  boundary  value 
problem  in  a  slab  [7] .  Results  of  the  analysis  of  the  slab  indicate  that 
under  high  rate  of  loading  the  coupling  term  in  the  energy  balance  equation 
should  not  be  ignored. 

This  paper  presents  the  numerical  method  and  some  results  obtained  in 
solving  the  coupled  dynamical  thermoelasticity  equations  in  a  long  hollow 
cylinder,  subjected  to  the  pressure  and  the  heat  flux  boundary  conditions 


at  the  inner  surface  and  ambient  environment  at  the  outer  surface  by  a  finite 
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element  method.  Several  numerical  schemes  are  studied  for  comparison  purpose. 
Besides  exhibiting  the  physics  of  the  problem,  this  paper  also  is  intended 
to  show  how  various  numerical  schemes  are  suited  for  these  problems. 


2.  Mathematical  Model 

The  classical  coupled  dynamical  thermoelasticity  equations  are  a  pair 
partial  differential  equations  governing  the  balance  of  the  linear  momentum 
and  the  energy  as  given  by  the  following  [2] . 

y72u  +.  (1  *  u)  grad  div  u  +  F  --  y  grad  T  *  ptt  (D 

and 

V2T  -  i  T  -  n  div  u  *  -  §•  (2) 

k  *  k 


where  u  is  the  displacement  vector  and  7  is  the  temperature,  X  and  y  are 
Lame's  constants,  p  is  the  density,  k  is  the  diffusivity,  Q  is  the  heat  source 
y  and  n  are  defined  as: 


y  ■  (3X  +  2y)a* 
n  -  YT0/pck 


(3) 


where  a*  is  the  coefficient  of  thermal  expansion,  c  is  the  heat  capacity,  Tq 
is  the  ambient  temperature. 


(5) 


Notice  that  due  to  polar  symmetry  the  displacement  vector  u  is  reduced  to  the 
radial  displacement  u,  whereas  the  tangential  component  is  identically  zero. 
Also,  it  is  assumed  that  F  =  0  ,  Q  =  0. 

In  Eqs.  (4)  and  (5) ,  we  have 

m  ■  y/(2 u  +  X)  and  Cj2«  (X  +  2u)/p  (6) 

The  initial  conditions  on  u  and  T  ares  ufr,0)-0,  TCr,0)-0. 

The  boundary  conditions  are: 

at  the  inner  boundary,  r  ■  ^ 


(X  +  2u)  ~  +  X  J  -  y(T  -  T0)  -  £Mt)  -  0 
and  (7) 

U  -  B*  [T  -  g*  (t)  ]  -  0 


at  the  outer  boundary,  r  ■  rQ 


(X  +  2u) 


X”  -  Y (T  -  T0) 


0 


and 


(8) 


£  +  32  T 


where  6*  *  hj/K  and  S2*  “  I12/K  ,  hj  and  hj  are  the  heat  transfer  coefficients 
at  the  inner  and  outer  boundaries  respectively.  K  is  the  thermal  conductivity 


Zt  is  noted  that  the  data  on  the  forcing  functions  f *  and  g*  are 
compatible  in  that  both  these  functions  are  continuous  at  the  initial  time. 


Introducing  the  following  set  of  non-dimensional  variables 


T  -  To 


To 


r . 

i 


-  u 

t  ■  and  u  »  — 

r-  r, 


(9) 


i  i 

and  after  some  algebra,  Egs.  (4)  and  (S)  ax re  transformed  into  the  following 
non-dimensional  equations  with  bars  removed. 


3  3u  3T 

lr(rlr)  +  7  +  +  X*  “ 


3r 


(10) 


and 


3  ,  3T 


-  +  crT  +  cl1j;(ru) 


3r 


(11) 


where  c  -  Cjr^/k,  XL  ■  kn/T0  and  X2  ■ 


mTr 


The  initial  conditions  are:  u(r,05-0,  T(r,0)-0  and  the  boundairy  conditions  are 
at  the  inner  boundary,  r  »  1 


3u 


~+a^-X2T-f(T) 


|“  -  BiT  +  8ig(T)  -  0 


(12) 


at  the  outer  boundary,  r  »  r 


3u  u  .  _ 
—  +  a  7  -  X2T 


(13) 


!£+  3*T 


where  0i  ■  Bi*^  ,  *  32*ri  »  °  *  X/(2u  ♦  X)  , 


g(T)  -  lg*(t)  -  T0]/To 


and  f(i)  ■  f* (t) / (X  ♦  2u) . 


3.  SPAC3  DISCRETIZATION 


Two  approaches  are  used  in  the  space  discretization  -  Galerkin  method  which 
uses  linear  shape  function  and  generates  consistent  mass  and  conductivity 
matrices,  and  central  explicit  method  which  uses  quadratic  shape  functions  and 
generates  lumped  mass  and  conductivity  matrices. 


3.1  Galerkin  Method.  In 
approximate  solutions  of  Eqs. 
expanded  in  terms  of  the  same 


this  method  it  is  assumed  that  both  of  the 
(10)  and  (11) ,  u  and  T  respectively,  can  be 
shape  function  ^(r)  (Figure  1). 


n 

u  ■  I  ip  .  (r)u.  (t) 
i-0 


-  n 

T  -  I  it.  (r)T  (T) 
i-0  X  1 


(14) 


where  n  is  the  number  of  nodes,  including  the  two  boundary  nodes,  u^  and 
are  the  approximate  nodal  values  of  the  displacement  and  the  temperature 
respectively. 

Substituting  Eq.  (14)  into  Eqs.  (10),  (11),  and  (12)  and  (13)  and  imposing 
the  condition 


(15) 


where  and  R2  are  the  residues  calculated  upon  substitution  as  stated,  and 
performing  integration  and  some  algebraic  work  lead  to  the  following  pair  of 
linear  algebraic  equations: 


[MU2] 0  +  [KU2] U  +  X2(KT2]I  -  F'  ' 


(16) 


and 

[CTllj  +  7  [KT1] I  +  XxCCUlig  -  F(1)  (17) 

C 

where  [MU2] ,  [KU2] , - ,  etc.  are  tridiagonal  N  x  N  matrices.  The  details 

of  deriving  the  solution  given  in  Eqs.  (15),  (16)  and  (17)  can  be  found  in 
Appendix  A. 


3.2  Central  Explicit  Method.  In  this  case  the  shape  functions  4k  (r) 
are  quadratic  instead  of  linear.  The  weighting  function  is  taken  to  be 
6 (r-r^) .  He  let  u  and  T  be  the  following  quadratic  expansion  within  each 
elesant. 


i+1 

u  -  I  iff.  (r)  u.  (T) 

i-1  1  1 

(18) 

i+1 

*  -  Z  <i>,  (r)  T  (T) 

i-1  1  x 


where  the  ij/'s  are  the  Lagrange  polynomials. 

Figure  2  illustrates  the  function. 

Using  the  same  notation  as  before,  we  impose  the  condition 

i+1 


Sfr-r^dr  ■  0  , 


i  -  0,1,2, - n 

j  *  1,2 


(19) 


The  respective  coefficient  matrices  of  Eqs.  (16)  and  (17)  for  this  procedure 
are  given  in  Appendix  B. 


4.  TIME  MARCHING 


Two  approaches  have  been  tried:  the  general  implicit  scheme  [8]  and  the 
three-point  recurrence  scheme  [9] • 

4.1  General  Implicit  Method  (GIM) .  For  convenience  we  shall  write  Eqs. 
(16)  and  (17)  in  a  single  equation 

MX  +  CX  +  KX  -  F  (20) 

where  X  -  {u,T}T 


AY  ■  3*  +  E 


(21) 


where 


V  *  ly  y  II 


F  -  {?(2)  9  F(1)}T 


MU  2 


1  0  1  0 


i  :  o 


0  CT1 


I 

J 


0  j-KU2  I  -X2KT2 


!  0 

I  “i 


-XjCUl)  0  ,  -  |  KT1 


Applying  GIM  to  Eq.  ( 21}  we  obtain  the  following  two-point  recurrence  scheme 


[A-&T0B] Y1+1» [A+At (1-0 ) B]Y  X+At [9?i+1+ (1-0  Jf1] 


(22) 


where  0  is  a  parameter,  0(3^1. 

The  selection  of  the  value  of  0  depends  on  the  problem  on  hand,  guided  by 
the  stability,  accuracy  and  economy  of  the  computation.  The  value  0  ”  0.667 
has  been  shown  to  be  a  good  choice  [9,10]. 


4.2  Three-Point  Recurrence  Method  (GTE) 

Let  the  shape  function  N^t)  be  defined  such  that 


-  -£  (1-0/2 

S,<T) 

-  (1-0  (l+£) 

23) 

■W1’ 

-  C(l+5)/2 

where 


£  -  t/At 
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Figure  3  illustrates  the  N^t)  functions. 


Let  X (t )  »  JJ.  .  (t)X.  , 
-  i-l  i-l 


N.  {T)X.  +  N .  .  (t ) X .  . 

1  1  1+1  1+1 


(24) 


The  residue  R  is  defined  as 


R*MX  +  CX  +  KX-F 


(25) 


We  require  that 


/ 


i+1 


R  W.  ( t ) dt  «  0 


(26) 


i-l 


where  w^t)  is  the  weighting  function.  Performing  the  operation  in  Eq.  (26)  yields 


[M  +  ydTC  +  SAt  K]Xi+1  +  (-2M  +  (1-2y)AtC 


(j  -  20  +  Y>AT2K]Xi  +  [M  -  (l-Y)ATC  +  (~  +  S-yJAt^JS^ 


-  FAt  «  0 


(27) 


where  8  and  y  are  two  parameters  depending  on  the  weighting  function  chosen. 


•A 


■MB  ^ld5 


Wi(5)d5 


-1 


and  1 

r 

y 


mj  Wi(C)(5+y)d Wi(C)dC 


■‘i 


r»-  y  ...... 
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'ect  of  various  combinations  of  3  and  y  for  the  solution  of  problems 
:junics  has  been  studied  [11]  . 

AND  DISCUSSION 

here  are  two  schemes  each  for  the  space  and  the  time  discretization, 
the  time  discretization  schemes  there  are  several  choices  of  values  of 
srs  9  ,  S  and  y>  many  combinations  of  the  space  and  time  discretizations 
e  for  the  numerical  solution  of  the  physical  problem. 

Dblems  were  selected  with  a  view  to  compare  the  results  of  computation 
g  the  relative  merit  of  each  scheme  regarding  the  stability  and 

sponses  to  a  Modified  Step  Stress  Input  (Wave  Equation  Only) . 
e  first  series  of  computation  was  done  to  test  the  wave  propagation 
When  the  temperature  T  was  removed  by  putting  X^  *  X2  *  0  and 
boundary  condition  deleted,  the  remaining  system  of  equations 
a  wave  equation,  for  which  analytical  solution  exists, 
te  computation  a  modified  step  stress  is  given  as 

f (T)  -  -0.1(1  -  e"10T) 

owing  computation  the  same  time  and  spatial  increments  are  used  for  all. 

ersus  T  Responses 

-  Figure  4  shows  the  result  of  computation  using  Galerkin  (GK)  in 
IM  in  time  discretization.  Three  different  values  of  9  were  used: 
ntral  difference),  9  -  0.667  (Gelerkin)  and  3*1  (backward).  The 
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figure  shows  the  radial  and  tangential  stresses  versus  time  at  the  midpoint 
(r  »  1.45).  The  results  indicate  that  in  the  case  of  S  *  0.667,  least  amount 
of  spurious  oscillation  exists,  whereas  the  other  two  schemes  of  GIM  lead  to 
either  spurious  oscillation  or  excessive  damping  as  evidenced  by  the  drastic 
reduction  in  amplitude. 

GK/GFS  -  Figure  5  is  the  same  set  of  responses  calculated  by  GK/GFS  with 
four  different  sets  of  values  for  3  and  y.  They  are  (1)  B  *  0.25,  y  *  .5 
(2)  S  *  1/6,  y  *  .5,  (3)  8  =  .5,  y  =  .6,  (4)  3  *  .8,  Y  *  1.5.  Results  show 
that  the  third  set  of  8  and  y  give  neither  spurious  oscillation  nor  excessive 
damping . 

CD/GIM  or  GFE  -  Figure  6  shows  the  computation  based  on  a  central  explicit  (CD) 
scheme  in  space  and  three  different  schemes  in  time  (1)  GIM,  9  »  0.667 
(2)  GFE,  6  *  .5,  y  -  .6  and  (3)  GFE,  6  -  .8#  Y  ■  1.0. 

Curve  (1)  shows  comparatively  little  dissipation  in  the  amplitude  responses, 
whereas  curve  (2)  and  (3)  exhibit  significant  dissipation. 

In  summary,  it  is  demonstrated  that  several  schemes  give  comparable 
results  in  the  wave  propagation  responses  which  are  accurate  and  stable, 
as  can  be  verified  at  least  qualitatively  by  a  method  such  as  characteristics. 

In  the  following,  the  fully  coupled  fields  in  the  cylinder  are  analyzed 
subject  to  either  stress  or  thermal  inputs  at  the  inner  bore.  In  the 
stress  input  we  use  both  pulse  or  step  type  input,  whereas  only  pulse  type 
input  is  used  for  the  temperature. 

In  this  paper  it  is  not  the  intention  to  give  details  regarding  the 
physical  problem  involved .  The  parameters  used  in  the  calculation 
derive  their  origin  from  the  application  problem.  For  reference  purposes, 
a  list  of  the  values  of  the  parameters  used  in  the  computation  is  given. 


The  amplitudes  of  the  pulses  and  the  step  are  realistic  values  corresponding 
to  the  physical  situation  in  the  application. 


-2t 

5.2  Responses  to  a  single  stress  pulse,  f(T)  »  -Q.015le  for  different 
r^/r^  -  At/Ay  -  and  6  -  values  (fully  coupled  equations) . 

The  next  series  of  three  figures  show  the  radial  stress  versus  time 
responses  of  the  cylinder  with  fully  coupled  equations  subject  to  a  single  stress 
pulse.  Figure  7  shows  for  a  cylinder  with  r^r^  ■  2.0.  Curve  (1)  ,  (3)  and  (4) 
correspond  to  GIM  with  different  values  of  8.  From  these  curves  it  appears 
that  when  9  *  0.667  as  shown  in  curve  (1)  the  result  is  the  best.  Curve  (3) 
with  6  ■  0.5  shows  too  much  oscillation,  and  curve  (4)  with  8  *  1.0  shows 
significant  "numerical  dissipation” . 

Curve  (1)  and  (2)  have  the  same  8-value  but  different  At/ Ay.  Curve  (1) 
with  At /Ay  «0.2  gives  better  result  in  that  it  has  less  "numerical  dissipation”. 

Figure  8  is  a  similar  computation  but  with  *  1.5 ,  8  ■  0.667,  but 

different  AT/Ay.  Curve  (3)  uses  a  very  small  At («  0.002),  thus  requires  too 
much  computation  time.  Curve  (2)  exhibits  significant  "numerical  dissipation". 
Curve  (1)  with  Ax/Ay  ■  0.2  is  the  best  compromise. 

Figure  9  is  calculation  for  a  very  thin  cylinder  with  yQ/y ^  *  1.1.  The 
schemes  use  8  *  0.667  and  Ar/Ay  -0.2  and  2.0  respectively.  For  AT/Ay  •  0.2 
the  result  shows  less  dissipation,  whereas  for  AT/Ay  ■  2.0  the  result  shows 
significant  "numerical  dissipation". 

Therefore,  from  the  above  results  we  chose  0  -  0.667,  At/Ay  »  0.2  for 
computing  responses  to  stress  inputs . 
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5.  3  Responses  to  a  Modified  Step  Stress,  £(x)»  -.003 (l-e~10T) 

Figure  10  shows  the  responses  in  the  radial  stress,  tangential  stress 
and  temperature  versus  time.  For  the  given  input  the  rise  time  is  0.8,  i.e. 
when  t*0.8,  f(T)*  -.003.  The  plot  is  for  the  midpoint  between  the  fifth  and  the 
sixth  nodes  with  the  non-dimensional  r=l. 45.  Since  the  non-dimensional  time 
t  is  scaled  against  the  travel  time  of  the  elastic  wave  through  a  distance  equal 
to  the  inner  radius  and  the  non-dimensional  distance  is  scaled  against  the  length 
equal  to  the  inner  radius,  it  would  take  a  unity  of  time  for  the  radial  wave  to 
travel  a  unit  distance  in  non-dimensional  scales.  Thus,  the  time  of  arrival 
of  the  radial  wave  should  be  about  r  »  0.45.  This  is  indeed  the  case  as  shown 
in  Fig.  10.  Notice  that  the  radial  wave  does  not  rise  to  the  maximum  value 
of  the  input  as  it  would  be  the  case  in  a  slab  [  7]  ,  due  to  the  effect  of  the 
curvature.  Observed  is  also  the  fact  that  the  reflected  wave,  being  tensile 
in  nature  due  to  the  free  surface  at  the  outer  boundary,  gives  rise  to  intervals 
of  time  when  the  radial  stress  drops  off  and  becomes  tensile  stress  for  a  short 
time.  A  periodicity  of  x-2  is  observed  from  the  figure  with  good  regularity 
within  the  time  plot.  The  travel  time  of  the  tangential  wave  should  be  about 
(2 it) (1.45)/1  or  9.11.  The  calculation  gives  a  good  approximate  check. 

The  temperature  response  in  the  cylinder  is  due  to  the  coupling  in  both 
field  equations.  Predictably ,  the  temperature  generated  is  small,  about 
0.6°C  *Vi  -1.4°C  referred  to  ambient  when  the  ambient  temperature  is  27®C.  It 
can  be  either  positive  or  negative  corresponding  to  a  rise  on  a  drop  from  the  ambient 
temperature  due  to  a  volumetric  compression  or  expansion.  This  effect  usually 
i3  very  small.  Here,  due  to  the  high  rate  of  change  of  strain,  the  coupling 


effect  is  not  ignorable. 


5.  4  Responses  to  a  Single  Triangular  Stress  Pulse 


A  triangular  pressure  pulse  of  duration  equal  to  2  and  amplitude  equal  to 
0.003  is  applied  at  the  inner  boundary.  Figure  11  shows  the  responses  vs 
time  and  Figure  12  the  distribution  of  responses. 

For  large  time  t  a  standing-wave  like  motion  is  observed  in  Figure  12 
as  each  point  goes  through  a  periodic  motion  of  different  amplitude ,  whereas 
the  two  boundaries  in  this  case  are  both  traction  free  after  x-2. 

After  x— 2 ,  the  maximum  non-dimensional  radial  stress  occurring  at 
x-2. 4  is  tensile  in  nature  and  equals  to  0.00194.  The  corresponding 
minimum  is  a  compressive  stress  .00164  at  T-7.4.  The  induced  temperature 
varies  from  1.2°C  at  t-7.4  to  -1.3*C  at  T-2.4,  referred  to  the  ambient 
temperature.  The  variation  is  about  ±5%  of  the  ambient  temperature .  These 
responses  are  not  graphed  in  the  figures  shown. 

5.5  Response  to  a  Temperature  Pulse,  g(x)  -  0.008155X0*°' 00051 

The  time  scale  for  the  temperature  response  is  about  four  decades  longer 

4 

than  the  stress  response.  The  responses  are  evaluated  at  times  when  x  ^  10  . 

In  the  numerical  computation  a  time  increment  was  chosen  with  the  same  factor 
in  mind. 

Figures  13  and  14  show  the  radial  and  tangential  stress  distribution 

4 

for  various  times.  Notice  that  when  x  is  10  the  physical  time  t  is  0.173  sec 
From  X  -  0.2  to  x  -  0.5  the  actual  time  elapsed  is  0.0519  sec.,  a  relatively 
short  time  for  heat  transfer . 


The  maximum  temperature  is  0.124  ('^64.2°C)  at  the  inner  surface  when 

4 

X  »  1.2  x  10  .  It  drops  down  quickly  as  shown  in  Figure  15. 


5 


6.  Summary 

In  summary,  several  observations  can  be  made  from  the  results  of  the 
numerical  experiments  of  this  investigation. 

(1)  it  is  possible  through  deliberate  experimentation  to  arrive  at  a 
feasible  numerical  scheme  to  limit  the  amount  of  numerical  dispersion 
and  dissipation  to  a  reasonable  level. 

(2)  The  coupling  introduced  in  the  energy  equation  could  cause  temperature 
fluctuation  due  to  a  rapidly  applied  stress  input  in  the  order  of  S 

to  10%  for  a  single  pulse  input,  depending  on  the  rate  of  application. 

(3)  The  effect  of  the  curvature  is  to  decrease  the  peak  response  of  the 
radial  stress,  in  the  presence  of  the  tangential  stress  component. 
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